library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.4     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(randomForest)
## randomForest 4.7-1.1
## Type rfNews() to see new features/changes/bug fixes.
## 
## Attaching package: 'randomForest'
## 
## The following object is masked from 'package:dplyr':
## 
##     combine
## 
## The following object is masked from 'package:ggplot2':
## 
##     margin
library(rfUtilities)
library(rfPermute)
## Welcome to rfPermute v2.5.2
## See rfPermuteTutorial() for a guide to the package.
library(ggpp)
## Registered S3 methods overwritten by 'ggpp':
##   method                  from   
##   heightDetails.titleGrob ggplot2
##   widthDetails.titleGrob  ggplot2
## 
## Attaching package: 'ggpp'
## 
## The following object is masked from 'package:ggplot2':
## 
##     annotate
library(dplyr)

#########
RF<-read.csv("/Users/jamesm./Desktop/THT_Data_5082.csv")
RF<-na.omit(RF)

RF=RF[,c(21,4:13)]


set.seed(123)
richness_rf <- randomForest(RF$Physical.Activity.from.Transportation..Score ~ ., data= RF,importance=TRUE,proximity=TRUE)
richness_rf
## 
## Call:
##  randomForest(formula = RF$Physical.Activity.from.Transportation..Score ~      ., data = RF, importance = TRUE, proximity = TRUE) 
##                Type of random forest: regression
##                      Number of trees: 500
## No. of variables tried at each split: 3
## 
##           Mean of squared residuals: 329.1773
##                     % Var explained: 59.36
richness_perm <- rf.significance(richness_rf, RF, nperm=99, ntree=501)
richness_perm
## Number of permutations:  99 
## p-value:  0.01 
## Model signifiant at p = 0.01 
##   Model R-square:  0.5858346 
##   Random R-square:  -0.1834711 
##   Random R-square variance:  0.01407312
set.seed(123)
richness_rfP<- rfPermute(RF$Physical.Activity.from.Transportation..Score ~ ., data = RF, ntree = 500,
                         na.action = na.omit, nrep = 99,num.cores = 3)

##
richness_dat <- importance(richness_rfP, sort.by = NULL, decreasing = TRUE)

#
richness_dat=as.data.frame(richness_dat)

##########
richness_dat$group2 <- row.names(richness_dat)

unique(richness_dat$group2)
##  [1] "Commute.Mode.Share...Auto..Score"       
##  [2] "Commute.Mode.Share...Auto..Raw.Value"   
##  [3] "Commute.Mode.Share...Transit..Raw.Value"
##  [4] "Commute.Mode.Share...Transit..Score"    
##  [5] "Commute.Mode.Share...Bicycle..Score"    
##  [6] "Commute.Mode.Share...Bicycle..Raw.Value"
##  [7] "Commute.Mode.Share...Walk..Score"       
##  [8] "Commute.Mode.Share...Walk..Raw.Value"   
##  [9] "Complete.Streets.Policies..Raw.Value"   
## [10] "Complete.Streets.Policies..Score"
#######################
richness_dat %>%
  as_tibble(rownames = "names") %>%
  data.frame() %>%
  mutate(label = if_else(X.IncMSE.pval < 0.001,"***",
                         if_else(X.IncMSE.pval <0.01,"**",
                                 if_else(X.IncMSE.pval<0.05,"*","  "))),
         X.IncMSE = as.numeric(X.IncMSE)) %>%
  mutate(group = if_else(label=="ns","In_sig","Sig"),
         names = forcats::fct_inorder(names)) %>%
  dplyr::arrange(X.IncMSE)  %>% 
  ggplot(aes(x = group2, y = X.IncMSE))+
  geom_bar(fill="tan2",stat = "identity",width=0.6,color="black")+
  geom_text(aes(y = X.IncMSE + 1,label = label),size=4)+
  labs(x = "", y = "%IncMSE", title="")+ 
  theme_classic()+ 
  #theme_bw()+
  coord_flip()+
  ####
  theme(
    panel.grid = element_blank(),
    panel.background =element_blank(),
    #legend.position = "none",
    legend.title = element_blank(),
    legend.text = element_text(color="black", 
                               size = 10, 
                               family = "serif",
                               face = "bold"),
    axis.title.y = element_text(size = 12, 
                                family = "serif",  
                                color = "black",
                                face = "bold"),
    axis.title.x = element_text(size = 12, 
                                family = "serif",  
                                color = "black",
                                face = "bold"),
    axis.text.x = element_text(size = 10, 
                               family = "serif",  
                               color = "black",
                               face = "bold"),
    axis.text.y = element_text(size = 10, 
                               family = "serif",  
                               color = "black",
                               face = "bold"))

RF<-read.csv("/Users/jamesm./Desktop/THT_Data_5082.csv")
RF<-na.omit(RF)

RF=RF[,c(21,16:19)]

set.seed(123)
richness_rf <- randomForest(RF$Physical.Activity.from.Transportation..Score ~ ., data= RF,importance=TRUE,proximity=TRUE)
richness_rf
## 
## Call:
##  randomForest(formula = RF$Physical.Activity.from.Transportation..Score ~      ., data = RF, importance = TRUE, proximity = TRUE) 
##                Type of random forest: regression
##                      Number of trees: 500
## No. of variables tried at each split: 1
## 
##           Mean of squared residuals: 194.4281
##                     % Var explained: 76
richness_perm <- rf.significance(richness_rf, RF, nperm=99, ntree=501)
richness_perm
## Number of permutations:  99 
## p-value:  0.01 
## Model signifiant at p = 0.01 
##   Model R-square:  0.7619681 
##   Random R-square:  -0.2068194 
##   Random R-square variance:  0.02194609
set.seed(123)
richness_rfP<- rfPermute(RF$Physical.Activity.from.Transportation..Score ~ ., data = RF, ntree = 500,
                         na.action = na.omit, nrep = 99,num.cores = 3)

##
richness_dat <- importance(richness_rfP, sort.by = NULL, decreasing = TRUE)

#
richness_dat=as.data.frame(richness_dat)

##########
richness_dat$group2 <- row.names(richness_dat)

unique(richness_dat$group2)
## [1] "Person.Miles.of.Travel.by.Walking..Raw.Value"        
## [2] "Person.Miles.of.Travel.by.Walking..Score"            
## [3] "Person.Miles.of.Travel.by.Private.Vehicle..Raw.Value"
## [4] "Person.Miles.of.Travel.by.Private.Vehicle..Score"
#######################
richness_dat %>%
  as_tibble(rownames = "names") %>%
  data.frame() %>%
  mutate(label = if_else(X.IncMSE.pval < 0.001,"***",
                         if_else(X.IncMSE.pval <0.01,"**",
                                 if_else(X.IncMSE.pval<0.05,"*","  "))),
         X.IncMSE = as.numeric(X.IncMSE)) %>%
  mutate(group = if_else(label=="ns","In_sig","Sig"),
         names = forcats::fct_inorder(names)) %>%
  dplyr::arrange(X.IncMSE)  %>% 
  ggplot(aes(x = group2, y = X.IncMSE))+
  geom_bar(fill="tan2",stat = "identity",width=0.6,color="black")+
  geom_text(aes(y = X.IncMSE + 1,label = label),size=4)+
  labs(x = "", y = "%IncMSE", title="")+ 
  theme_classic()+ 
  #theme_bw()+
  coord_flip()+
  ####
  theme(
    panel.grid = element_blank(),
    panel.background =element_blank(),
    #legend.position = "none",
    legend.title = element_blank(),
    legend.text = element_text(color="black", 
                               size = 10, 
                               family = "serif",
                               face = "bold"),
    axis.title.y = element_text(size = 12, 
                                family = "serif",  
                                color = "black",
                                face = "bold"),
    axis.title.x = element_text(size = 12, 
                                family = "serif",  
                                color = "black",
                                face = "bold"),
    axis.text.x = element_text(size = 10, 
                               family = "serif",  
                               color = "black",
                               face = "bold"),
    axis.text.y = element_text(size = 10, 
                               family = "serif",  
                               color = "black",
                               face = "bold"))

###################

us_states <- map_data("state")
election<-read.csv("/Users/jamesm./Desktop/THT_Data_5082.csv")
us_states_elec <- left_join(us_states, election,by="region")

library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
#
p0 <- ggplot(data = us_states_elec,
             mapping = aes(x = long, y = lat,
                           group = group, fill = us_states_elec$Physical.Activity.from.Transportation..Score))+
  geom_polygon(color = "gray90", size = 0.1) +
  coord_map(projection = "albers", lat0 = 39, lat1 = 45) +
  theme(legend.title = element_blank())+
  scale_fill_gradientn(colours = RColorBrewer::brewer.pal(11, "RdBu"))+
  labs(title = "Physical Activity from Transportation Score",fill="")+
  theme_bw()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
ggplotly(p0)
## Warning: Use of `us_states_elec$Physical.Activity.from.Transportation..Score` is
## discouraged.
## ℹ Use `Physical.Activity.from.Transportation..Score` instead.
###
p0 <- ggplot(data = us_states_elec,
             mapping = aes(x = long, y = lat,
                           group = group, fill = Commute.Mode.Share...Auto..Raw.Value))+
  geom_polygon(color = "gray90", size = 0.1) +
  coord_map(projection = "albers", lat0 = 39, lat1 = 45) +
  theme(legend.title = element_blank())+
  scale_fill_gradientn(colours = RColorBrewer::brewer.pal(11, "RdBu"))+
  labs(title = "Commute Mode.Share Auto Raw Value",fill="")+
  theme_bw()



ggplotly(p0)
####
p0 <- ggplot(data = us_states_elec,
             mapping = aes(x = long, y = lat,
                           group = group, fill = Commute.Mode.Share...Auto..Score))+
  geom_polygon(color = "gray90", size = 0.1) +
  coord_map(projection = "albers", lat0 = 39, lat1 = 45) +
  theme(legend.title = element_blank())+
  scale_fill_gradientn(colours = RColorBrewer::brewer.pal(11, "RdBu"))+
  labs(title = "Commute Mode Share Auto Score",fill="")+
  theme_bw()


ggplotly(p0)